Stacked clusters of polycyclic aromatic hydrocarbon molecules 
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Clusters of polycyclic aromatic hydrocarbon (PAH) molecules are modelled using explicit all-atom 
potentials using a rigid body approximation. The PAH's considered range from pyrene (CioHg) to 
circumcoronene (C54H18), and clusters containing between 2 and 32 molecules are investigated. In 
addition to the usual repulsion-dispersion interactions, electrostatic point-charge interactions are 
incorporated, as obtained from density functional theory calculations. The general electrostatic 
distribution in neutral or singly charged PAH's is reproduced well using a fluctuating charges anal- 
ysis, which provides an adequate description of the multipolar distribution. Global optimization is 
performed using a variety of methods, including basin-hopping and parallel tempering Monte Carlo. 
We find evidence that stacking the PAH molecules generally yields the most stable motif. A struc- 
tural transition between one-dimensional stacks and three-dimensional shapes built from mutiple 
stacks is observed at larger sizes, and the threshold for this transition increases with the size of 
the monomer. Larger aggregates seem to evolve toward the packing observed for benzene in bulk. 
Difficulties met in optimizing these clusters are analysed in terms of the strong anisotropy of the 
molecules. We also discuss segregation in heterogeneous clusters and vibrational properties in the 
context of astrophysical observations. 



I. INTRODUCTION 



Polycyclic aromatic hydrocarbons (PAH's) have been 
proposed as the carriers of a family of interstellar aro- 
matic infrared bands (AIB's), observed in many astro- 
nomical objectSfii^ The most intense of these bands, 
found near 3030, 1610, 1315-1282, 1150 and 885 cm"! 
respectively, are observed systematically from different 
regions of the interstellar medium heated by starlight. 
The intensity of these features suggests that PAH's are 
the most abundant complex polyatomic molecules in the 
interstellar medium, accounting perhaps for as much 
as 20% of all carbon in our galaxy,'^ Many laboratory 

studiea'^i^i^i''i^i^i-^°i-^-^i-^^i-^'^i-^'^ and quantum chemical calcu- 
lations or modelal^ll'^ll7ll8,19l20,21,22^23,24 of gi^glg PAH 

molecules have been performed in the past, but none 
has yet provided convincing evidence that single PAH 
molecules are actually present in the interstellar medium. 
This result is likely due to the specific nature of the in- 
terstellar species, in relation to their formation mech- 
anism and further processes such as photodissociation. 
Boulanger et alm^ and Bernard et alm^ suggested that 
the observed free-flying PAH's are produced by photoe- 
vaporation of larger grains. Cesarsky et al^^ attributed 
the presence of an infrared continuum observed in the 
reflection nebula CED 201 to very small carbonaceous 



grains eventually leading to the AIB's carriers. In recent 
work, Rapacioli et alm^ presented strong evidence that 
these grains are PAH clusters, and a typical lower size of 
400 carbons per cluster was inferred. 

Generally speaking, very little is known about 
the structure of PAH clusters. In comparison, as- 
semblies of benzene molecules have received much 
more attentioni^9i30,3ii32,33,34,35,36,37,38,39,40,4i At empir- 
ical levels of theory, van de WaalSS*^ and the groups 
of Stacc,3i Bartell2& and Whetten and Eastev^iL^^ 
have investigated assemblies containing up to 13 ben- 
zene molecules. These studies showed a marked prefer- 
ence for the Wefelmeier growth scheme^ also observed 
in argon clusters."^ The (benzene) 13 cluster, in particu- 
lar, was shown to exhibit very stable icosahedral-based 
structures some of which have received experimental 
support. At more sophisticated levels, electronic struc- 
ture calculations have also been performed by several 
authors .Si^iiS42i^ The cationic species (C6H6)+ have 
been investigated experimentally from ion mobility mea- 
surements by Rusyniak et almi- These authors also per- 
formed structural optimization using the OPLS force 
field.i^ 

Naphthalene clusters have been studied by van de 
Waal,^^ who identified patterns similar to benzene clus- 
ters. More recently, the structure of the naphtha- 
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lene trimer obtained from experimental studies^ was 
shown to agree well with one derived from ab initio 
calculations!^ 

Anthracene clusters containing up to 5 molecules were 
investigated theoretically and experimentally by Piuzzi 
and coworkers4£ Despite their relatively small size, these 
clusters seem to show quite different structures, forming 
a mainly two-dimensional pattern in which the long axes 
of the molecules are aligned, their centers of mass lying in 
the same perpendicular plane. Song and coworkers^- re- 
ported mass spectrometry measurements for anionic clus- 
ters containing up to 16 anthracene molecules. 

Heterogeneous clusters made of one or several benzene 
molecules and a single naphthalene, anthracene, pery- 
lene or tetracene molecule have been studied, sometimes 
in the presence of a solventi^fii^ Data for larger PAH 
molecules is significantly more scarce. The coronene 
and circumcoronene dimers were investigated by Miller 
and coworkers. Aggregates of PAH molecules were as- 
sumed by Seahra and Dulej*S to form stacks. These 
authors used data from graphite to estimate some trans- 
lational vibrational modes of stacked PAH's. However, 
this assumption did not employ any atomistic modelling. 
MarzeC)^ using an extension of the MM2 force field^^ 
and semi-empirical quantum mechanical approaches, per- 
formed local optimizations of stacked structures contain- 
ing up to 6 molecules. In this study, the stacks were 
shown to be stable for PAH's ranging from dibenzopy- 
rene to dicoronene.^'* Perlstein— investigated the transi- 
tion from one-dimensional stacks to monolayers and crys- 
tal packings of various molecular species, also using the 
MM2 classical potential. 

While bulk polyaromatic compounds are very difficult 
to study experimentally, nearly pure macroscopic sam- 
ples of hexabenzocoronene molecules were shown to form 
a twinned crystal. Crystallization via epitaxial growth 
over a graphite surface has also been reportedi^ Com- 
puter simulations results by Khanna and coworkersiSi. 
indicate that large sets of coronene or circumcoronene 
molecules do indeed crystallise and melt through a first- 
order process without exhibiting intermediate liquid crys- 
tal behavior. More generally, the bulk phases of PAH 
molecules can be expected to show similar phase dia- 
grams as the model discotic particles studied by several 
authorsi^SiiS^ 

Our present interest lies in the size range between the 
very small and very large sizes, for aggregates contain- 
ing about 1000 carbon atoms. We wish to address the 
following questions: 

(i) Are the stacks of PAH molecules actually the most 
stable form for clusters? 

(ii) How large can these stacks be as a function of the 
PAH monomer size? 

(iii) How does the cluster structure evolve toward the 
bulk morphology? 



In addition to these main concerns, we intend to address 
not only clusters of the same PAH molecule, but also 
some heterogeneous assemblies, as the experimental sit- 
uation we are referring to does not have restrictions on 
the variety of PAH's involved. 

Locating the most stable structures of molecular clus- 
ters can be quite difficult i^^iS^*^ For example, water clus- 
ters have been found to exhibit a multi-funnel potential 
energy surface ;^ '^'^^'^^ where numerous low energy min- 
ima exist, but can only be interconnected by long path- 
ways. Global optimization is significantly harder for rigid 
body water clusters than for most atomic clusters with a 
comparable number of degrees of freedom. ^■^ PAH clus- 
ters, which are made of very anisotropic monomers, rep- 
resent a novel challenge for molecular optimization, some- 
where in between the difficulties associated with atomic 
clusters and those of biological molecules. 

In the present work, we have attempted to locate the 
stable structures of assemblies of relatively large PAH 
molecules ranging from pyrene to circumcoronene. We 
combine results from electronic structure calculations 
and atomistic modelling to describe the intermolecular 
interactions, within a rigid body approximation for the 
PAH molecules. Global optimization using the basin- 
hoppingSiZS and parallel temperingSi. Monte Carlo meth- 
ods then provide estimates of the most stable isomers, for 
clusters containing up to 32 molecules. 

The paper is organized as follows. In the next section, 
we give the basic description of our model, as well as some 
details about the electronic structure calculations and the 
optimization procedure. Section IlIII presents our results 
for the structures, starting with the dimers. A more de- 
tailed analysis is provided in the case of coronene clusters, 
and the relative stability of some important structural 
motifs is investigated. Sec. Illll ends with some results for 
heterogeneous clusters. In Sec. II VI we briefly discuss in- 
termolecular vibrational modes. We comment in Section 
0on the difficulties encountered in our global optimiza- 
tion approach and in the robustness of our results with 
respect to the quality of the atomistic modelling. Finally, 
we summarize and discuss the astrophysical relevance of 
our results in Sec. IVII 



II. METHODS 

Our goal here is to provide good candidates for the 
most stable structures of clusters of PAH molecules. 
These structures are expected to be relevant at very low 
temperatures, where intramolecular vibrations are un- 
likely to be excited. This legitimises our main approxi- 
mation that the PAH molecules can be treated as rigid 
bodies. 
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A. Intermolecular potentials 

The PAH molecules considered are pyrene (CieHio), 
coronene (C24H12), ovalene (CaoHig), hexabenzo- 
coronene (HBC, C42H18), octabenzocoronene (OBC, 
C46H18), and circumcoronene (C54H18). The cluster sizes 
range from 2 to 32 depending on the monomer, so we will 
be dealing with a rather large number of atoms. Con- 
sidering the known difficulty of global optimization's^ it 
seems obvious that the atomistic modelling should be 
limited to a moderate numerical cost — but retaining the 
important chemical features. In particular, the multipo- 
lar electrostatic description employed in Ref. is not ap- 
propriate for an initial survey, although such ideas could 
provide useful corrections. Following previous efforts on 
benzene clusters by several groups we have 
chosen to describe the intermolecular potential in PAH 
clusters as the result of two contributions: 



i<j a£i /3ej 



) + ^Q(r^„ (1) 



where Vlj (r^^ ) denotes the repulsion-dispersion energy 
between atom of molecule i and jp of molecule j, 
and Vq the electrostatic interaction between the partial 
charges carried by the molecules. Here and in the fol- 
lowing, the PAH molecules are denoted by roman letters, 
and greek letters are used for the atoms of each molecule. 

The Lennard- Jones (LJ) form is used for 
the repulsion-dispersion, with parameters 

(ctcc,o-hh,o-ch,£cc,£hh,£ch) taken from Ref. |2^ 
The electrostatic contributions are based on the point 
charges {qi} located on the atom sites. 





s 

Exp. Cont. 


P 

Exp. Cont. 


d 

Exp. Cont. 


c 


2.382013 -0.242140 
1.443065 0.185265 
0.405847 0.591283 
0.138427 1.0 


8.609570 0.043653 
1.943550 0.209497 
0.542798 0.502761 
0.152496 1.0 


0.80000 1.0 


H 


13.24876 0.019255 
2.003127 0.134420 
0.455867 0.469565 
0.124795 1.0 


0.90000 1.0 





TABLE I: Atomic GTO basis sets employed for carbon and 
hydrogen. 



faces the overlap problem and the corresponding results 
may depend strongly on the basis set. Other definitions 
of atomic charges can be found, such as the natural bond 
orbital (NBO) definitionZ^ or the Bader definition^ in 
which the charges are determined from the topology as- 
sociated with the density. For the present purpose how- 
ever, since the charges are essentially used to extract an 
intermolecular potential rather than describing the in- 
tramolecular bonding, we choose to compute the charges 
which best fits the actual electrostatic potential on the 
van der Waals surface associated with each molecule. 
This choice is similar to the strategy followed by Piuzzi 
and coworkers,*^ which provided reasonable results for 
the anthracene dimer by consistently reproducing the ab 
initio interaction energies for various relative orientations 
of the molecules. 



B. DFT calculations and electrostatic modelling 

The molecular geometries and partial charges were ini- 
tially determined from electronic structure calculations 
on PAH monomers. The molecular geometries were those 
determined in tight-binding calculations. The partial 
atomic charges were initially determined from electronic 
structure calculations using the B3LYP hybrid functional 
and the GaussianOS packagen^S Those calculations were 
carried out using core pseudopotentials on carbon as de- 
termined by Durand and Barthelat^ and tabulated in 
semi-local form by Bouteiller et ali^ The gaussian type 
orbital (GTO) basis set used was [4s/4p/ld] contracted 
into [31/31/1] on the carbon atoms and [4s/lp] con- 
tracted into [31/1] on the hydrogen atoms. The basis 
sets are specified fully in tabled 

In addition to the neutral PAH molecules, calculations 
were also performed on the singly charged anions and 
cations at the neutral geometry. These results allow us 
to investigate the relative stability of zwitterionic species. 

The extraction of point atomic charges from ah initio 
calculations is rather problematic since charges are not 
uniquely defined. The well-known MuUiken definition^^ 



The fitted charges were determined here using the 
GaussianOS package for all the species mentioned above. 
As expected, the electrostatic potential fitted (EPF) 
charges are significantly different from those given by the 
MuUiken analysis. They show more continuity from one 
PAH to one another, especially in terms of the popula- 
tions on hydrogens with respect to carbon atoms. With 
the present basis set, the MuUiken analysis often leads to 
negative charges on hydrogen atoms, whereas both NBO 
and EPF charges for neutral PAH's are always positive 
on hydrogen. Table HTl illustrates the values obtained for 
coronene. The other PAH's behave similarly, except that 
they show larger fiuctuations within a given molecule, 
especially when the symmetry is lower. 

After the charges were determined, a simple electro- 
static model was constructed to rationalize the DFT re- 
sults. The model is similar to the charge equilibration 
scheme of Mortier et al.^ and also to the fluctuating 
charges model of Rappe and Goddardi2^ Interestingly, 
some links between DFT and fluctuating charges poten- 
tials have been established]^ Briefly, at a given molecular 
geometry the atomic charges {qi\ are defined in order to 
minimize the global electrostatic energy under the con- 
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shell 


MuUiken 


NBO 


EPF 


Cl(6) 


0.0040 


-0.0082 


0.0004 


C2(6) 


-0.0295 


-0.0492 


0.0799 


C3(12) 


0.0299 


-0.1925 


-0.1697 


H(12) 


-0.0172 


0.2212 


0.1295 



TABLE II: MuUiken, natural bond orbital, and electrostatic 
potential fitted charges for carbon and hydrogen atoms of 
coronene. The carbon atoms are labelled according to their 
radial location from central to peripheral, and the number of 
equivalent atoms is also indicated in parentheses. 

straint of total charge conservation: 

i ^ i<j \ i / 

(2) 

In the above equation, and Ui are the electronega- 
tivity and hardness of atom i, respectively, and Jy = 
(^fj +7ij^)~^^^- The Lagrange multiplier A ensures that 
the total charge is constant and equal to Qq. This elec- 
trostatic model has only 4 parameters, namely ec — £Hj 
Uc = 7CC7 = 7hh, and 7ch- These parameters were 
optimized to reproduce the partial charges of the neutral 
coronene molecule, as obtained from density functional 
calculations. We found the charges carried by the atoms 
of all other PAH molecules to be very satisfactorily re- 
produced by this model, always within 5%, not only for 
neutrals, but also for anionic or cationic PAH's. There- 
fore, the complete electrostatic atomistic description is 
gathered into the parameters of this simple model, which 
read £c - eh = 1-86 eV, Uq = 0.98 eV, Uh = 1.65 eV, 
and 7CH = 2.11 eV. 

Before going on, it is important to notice that the 
present fluctuating charges model not only accounts for 
(first-order) Coulombic interactions. The charges ob- 
tained are intrinsically the solution of a linear many-body 
equation, and the resulting energy incorporates polariza- 
tion effects. In addition to providing all the important 
electrostatic characteristics of the single molecule, the 
fluctuating charges model can be extended to treat an 
assembly of PAH molecules. For this purpose, Eq. Q 
must be extended to include interactions between all the 
atoms, and the conservation of total charge over each 
molecule is achieved through the introduction of a cor- 
responding number of Lagrange multipliers^ In princi- 
ple, charge transfer among a limited or complete set of 
molecules could also be investigated, but we have not 
considered this possibility here. 

Allowing charges to fluctuate within a PAH cluster has 
a very heavy computational cost, which turned out to 
be incompatible with performing large scale optimiza- 
tion. We therefore restricted the model to flxed charges. 
However, for some of the typical equilibrium structures 
we found, we checked that the assembly of several PAH 
molecules was consistent with the flxed charges approxi- 



mation. To do this we computed the average {\qia — Qia \) 
over all atoms and all molecules from the charges {q^^ } in 
the single molecules and their values {qia} in the relaxed 
cluster. For all clusters investigated, the average charge 
always remained lower than 0.1%. Obviously, at flnite 
temperatures intramolecular vibrations are likely to in- 
duce signiflcant variations in the effective atomic charges, 
and the flxed charges approximation will become less ac- 
curate. 



C. Global optimization 

We attempted to locate good candidates for the global 
minima of PAH clusters using several methods. The 
basin-hopping (BH), or Monte Carlo plus minimization 
methodSSiZS, was used, especially for the smaller clusters. 
This unbiased algorithm performs a random walk on the 
transformed potential energy surface, which consists only 
of the local minima. This walk requires quenches to be 
carried out after each molecular displacement. Here, a 
displacement means that all molecules are rotated and 
translated randomly by some (rather large) amount, at 
the same time. 

The basin-hopping method was previously used for a 
variety of other atomic and moleculariSiSiSi clusters. 
However, we found it rather inefflcient for moderately 
large PAH clusters, for reasons that will be discussed 
below. Therefore, we also employed parallel tempering 
Monte CarloSi as an alternative optimization tool. Even 
though this method was originally designed for reduc- 
ing broken ergodicity in glassy systems, it has since been 
used successfully in global optimization problems It 
essentially consists in performing several MC trajecto- 
ries at increasing temperatures, and allowing exchange 
moves between adjacent trajectories in addition to the 
standard molecular displacements. In practice, we chose 
20 temperatures in the range 1 K< T <80 K. Periodic 
quenches from some trajectories were also considered to 
provide extra isomers. Finally, many low energy struc- 
tures were found by construction from simple motifs, us- 
ing the structures obtained from unbiased methods for 
guidance. 



III. STRUCTURES 

The interaction between molecules is essentially ad- 
ditive within our model. However, due to the strong 
aspherical character of the molecules, this interaction is 
also expected to be highly anisotropic. In view of the 
numerous previous results for dimers of PAH's, we flrst 
focus on these species as a flrst step toward larger aggre- 
gates. 
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A. Dimers 

The coronene molecule provides a convenient ex- 
ample of large PAH's, especially from the computational 
point of view. The most stable structures of its dimer are 
shown in Fig. ^ (a) and (b). These structures are per- 




(a) E = -94.90 kJ/mol (b) E = -93.66 kJ/mol 




(c) E^ = -92.35 kJ/mol 

FIG. 1: Low energy structures of the coronene dimer. (a) 
Twisted stack; (b) parallel-displaced stack; (c) superimposed 
stack (saddle point). 

feet stacks, in the sense that the sixfold axes and molec- 
ular planes are parallel. A detailed analysis of the var- 
ious contributions to the binding energy shows that the 
repulsion-dispersion between carbon atoms dominates. 
Therefore it may be energetically favorable to shift the 
molecules along some parallel displacement, or to rotate 
them around their common sixfold axis. The resulting 
structures show an increase in the number of nearest 
neighbors, but the parallel displaced geometry is slightly 
penalized at its edges due to fewer van der Waals bonds 
involving hydrogen atoms. Interestingly, the perfectly 
superimposed stack, Fig^, is not a true minimum but 
a saddle point for our potential. This result contrasts 
with that obtained by MarzeC)^ who found the face-to- 
face configuration to be the lowest in energy. It should be 
noted here that the charges used by this author, obtained 
from the MIND03 and ZINDOl semi-empirical quantum 
mechanical calculations, are significantly lower than ours 
(always between 10"'^ and 10^^). However, even after 
removing the electrostatic contribution from our calcula- 
tion we still find that this conformation is a saddle. This 
result may indicate that the sirnple steepest-descent local 
minimization performed in Ref. 54 did not fully converge. 

A pair of pyrene molecules also exhibits stacks, as 
shown in Fig. |21 The superimposed structure analogous 
to (c) of Fig. n] spontaneously transforms into a twisted 
stack where the angle between the long axes is about 
30°. Another stable stack is also found with perpendicu- 
lar long axes. This isomer, (b) in Fig.|21 will be denoted 




(a) E = -54.47 kJ/mol (b) E = -53.24 kJ/mol 

FIG. 2: Low energy structures of the pyrene dimer. The angle 
between the long axes is (a) 50° and (b) 90°. 



as 'staggered' below. 

Smaller aromatic molecules also display stable stacked 
dimers, even though T-shaped geometries are found com- 
petitive for benzene, naphthalene and anthracene4^*^2iS 
(Benzene)2 is slightly more stable when T-shaped, but 
increasing the number of aromatic rings favors stacked 
shapeSf^- In both (naphthalene)2 and (anthracene)2, the 
crossed stacks are marginally preferred to the parallel 
displaced stacks^^i^ However, we note that, among T- 
shaped conformations, those with parallel long axes show 
an enhanced stability42i 

The trends seen in the naphthalene and anthracene 
dimers are also observed in the present work for 
(coronene)2 and (pyrene)2. In the latter case, the paral- 
lel displaced isomer leads to structure (a) in Fig. Rafter 
minimization. In the circumcoronene dimer, not shown, 
the crossed stacked structure is also the most stable, the 
parallel displaced stack being 1.7 kJ/mol higher in en- 
ergy. Again, the perfectly superimposed dimer is not a 
stable minimum. 

The favored dimer structure results from a balance be- 
tween maximizing the number of weak C-C intermolecu- 
lar contacts, and minimizing the Coulomb repulsion be- 
tween hydrogen atoms. The gain in repulsion-dispersion 
energy is larger than the Coulomb penalty in the crossed 
stack relative to the parallel displaced stack. We have 
also investigated cohesion in zwitterionic (PAH+PAH") 
stacks. For all the clusters investigated, the crossed 
stacks were further stabilized with respect to the parallel- 
displaced isomer. In these cases, the Coulomb interaction 
becomes strongly attractive instead of being slightly re- 
pulsive for the neutral forms. 

/,From the DFT calculated vertical ionisation poten- 
tials and electron affinities, and the binding energy of 
the ion pair, the energy of the zwitterion can be es- 
timated relative to the neutrals. The binding energy 
was calculated using the same intermolecular potential, 
with fixed partial charges corresponding to each partic- 
ular ion. This scheme neglects the resonance interaction 
PAH+PAH- ^PAH-PAH+. We find the relative zwit- 
terion energies AE = 4.23, 3.38, and 2.87 eV above the 
the neutral minimum for the pyrene, coronene, and OBC 
dimers, respectively. These values are similar to the re- 
sult reported by Piuzzi et al. for anthracene trimers of 
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around 3.5 eVi*^ 

Finally, the distance between molecular planes was 
found to lie between 3.56 A for (pyrene)2 and 3.50 A for 
(circumcoronene)2. These values compare well with the 
data reported by Marzecp'?^ as well as the experimental 
value for crystal packings obtained by Goddard et al. on 
HBC,^ 



B. Coronene clusters 

The lowest-energy structures of coronene aggregates 
are found to be single stacks as long as the number of 
molecules does not exceed eight. In these stacks, the ori- 
entations alternate so that adjacent molecules are stag- 
gered. Single stacks remain local minima for larger num- 
bers of molecules, mainly due to the short-range charac- 
ter of the interaction. However, above eight molecules 
lower energy structures appear, as illustrated in Fig. O 
for (coronene) 10- 




(a) (b) (c) 

-926.42 kj/mol -893.63 kj/mol -882.41 kj/mol 



FIG. 3: Lowest energy structures of the (coronene) lo cluster. 

The most stable conformation is shown in Fig. In 
this global minimum the two 'C '-shaped stacks adopt a 
configuration that looks rather like a handshake. The 
tilted stacks of Fig. |3J:, which we expected to be natural 
evolution of the one-dimensional stack, turn out to be 
less favorable at this size. The special stability of the 
'handshake' structure is due to its compactness, while 
conserving the primary stack motifs as basic units. Not 
surprisingly, other combinations of stacks such as 4-f6 
or 3-f 7 are higher in energy than the 5+5 structure of 
Fig.Et. 

Building upon the (coronene) lo cluster, extra 
molecules first tend to add at the tips of existing stacks. 
However, the 'handshake' 8-1-8 stacks are actually less 
stable than a 8-1-4-1-4 triple stack where the two short 4- 
stacks make opposite turns around an S-shaped 8-stack. 



In comparison, the 6-1-5-1-5 triple stack built by adding 
a third column to the structure of Fig. O; in a triangu- 
lar fashion is slightly higher in energy. This difference is 
mainly due to the geometrical frustration arising when 
trying to alternate the new molecular planes with those 
already present. 

We attempted to grow the S-shaped structure fur- 
ther by adding molecules at the six extremities, but this 
proved not to be the optimal route for larger clusters. 
There are many ways to put together 4-stacks in order 
to generate a (PAH)32 cluster. The four most stable such 
conformations are sketched in Fig. 0] 




(a) -3514 (b) -3490 




(c) -3487 (d) -3474 

FIG. 4: Lowest energy stacked structures of the (coronene)32 
cluster. The total binding energies are given in kj/mol. 

Among these conformations, three schemes based on 
alternating stacks oriented perpendicularly are very close 
in energy, independent of the stack sizes. However, four 
8-stack, tilted columns with parallel axes provide the 
most stable conformation at this size. This structure, 
depicted in Fig. appears as a clear precursor to the 
herringbone crystal packing reported by Khanna et al^ 
from computer experiments. Even though we did not 
carry out a systematic search for minima above 20 PAH 
molecules, decreasing the stack size from the (coronene)32 
herringbone motif of Fig. 0Ji remained the most stable 
structure down to less than 24 molecules. Below this 
size, perpendicular stacks similar to Fig. ^ become the 
global minimum. The fully intertwinned structure ex- 
emplified in Fig. ^ is only optimal at sizes lower than 
20. Therefore, it seems that bulk structure appears in 
coronene at about ~ 24 molecules. 

The appearance of bulk features in the structural prop- 
erties of clusters has attracted much attention in the past, 
especially in molecular clusters JS^-^^'®" It was noticed in 
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particular that more isotropic species such as argon or 
nitrogen exhibit the bulk crystalline character at much 
larger sizes than anisotropic molecules like C02i^ The 
results for coroncne thus seem to confirm this observa- 
tion. 



C. Stability 

The previous results obtained for coronene clusters 
generally hold for other PAH molecular units. The rel- 
ative stability of single- or multiple-stack structures is 
shown in Fig.|3 Pyrene clusters form single stacks, then 
'handshake' stacks above seven molecules. In circum- 
coronene clusters, the one-dimensional stack remains the 
global minimum until 17 molecules are reached. The sin- 
gle stacks remain stable longer for larger PAH's, prob- 
ably because this motif is favorable until the length ex- 
ceeds the diameter of the molecule itself. From = 18 
and above, circumcoronene clusters favor the herringbone 
conformation pictured in Fig. |2t for coronene. This re- 
sult again confirms the idea that increasingly anisotropic 
molecules show earlier signatures of bulk structure. 
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The relative stability of stacked configurations is also 
addressed in Fig. where we compare the dissociation 
energy of double stack structures (herringbone or 'hand- 
shake' conformations) into their two stacks with that of 
the most weakly bound molecule. The latter molecule 
is located at one of the four extremities of the stacks. 
Again, we consider pyrene, coronene and circumcoronene 
clusters. Obviously, the dissociation energy of a double 
stack into two stacks grows approximately linearly with 
stack size. Conversely, the binding energy of one single 
molecule remains nearly constant, and depends mainly 
on the monomer. These trends, evident in Fig. El ex- 
plain why it becomes more interesting to remove a sin- 
gle molecule only after some crossover size, below which 
the extra stability of the small stacks dominates. The 
crossover sizes, seven for pyrene, nine for coronene and 
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17 for circumcoronene, correlate with the change in sta- 
bility of double stacks with respect to single stacks. This 
result emphasizes that, unlike most weakly bound clus- 
ters, monomer dissociation may not always be the lowest 
energy dissociation pathway for such species. 



D. Demixing in heterogeneous clusters 

To date, no specific PAH has been clearly identified 
as the carrier of the interstellar AIB's. In fact, these 
bands probably result from the emission of a popula- 
tion of several types of PAH molecules jiSiiS and we ex- 
pect interstellar PAH assemblies to be quite heteroge- 
neous. It is therefore important to characterize the ef- 
fects of mixing on the most stable conformations. Ta- 
ble mil shows the relative stability of the various pos- 
sible stacks containing two pairs of identical molecules 
based on sixfold symmetry. As in homogeneous clus- 
ters, we were not able to locate structures more stable 
than these perfect stacks. Stack (a), in which the larger 
PAH's occupy interior sites, maximizes the number of 
nearest neighbors between carbon atoms and is ener- 
getically the most favorable. Similar arguments explain 
the ordering of isomers (a-d). As expected, the relative 
stability increases with the disparity between the PAH 
molecules: the energy of stack (b) is smaller than those of 
stack (a) by about 10% for (coronene-fcircumcoronene)2, 
but only 8.7% for (coronene-|-HBC)2 and 3.4% for 
(HBC+circumcoronene)2. 

This segregation property also appears for larger clus- 
ters. As an example we have investigated the hetero- 
geneous assemblies containing 6 large and 4 small PAH 
molecules of the coronene, HBC, and circumcoronene 
types. The putative global minima we found are all based 
on the same structures obtained for the homogeneous 
coronene cluster depicted in Fig. O For each of these 
isomers, the optimal locations of the smaller PAH are at 
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stack type 










(a) 


(b) 


(c) 


(d) 


(coronene)2 (HBC)2 


-473.3 


-435.5 


-515.5 


-376.4 


(coronene) 2 (circ. ) 2 


-572.8 


-518.7 


-478.0 


-418.4 


(HBC)2(circ.)2 


-722.3 


-695.9 


-687.4 


-663.1 



TABLE III: Binding energies (kj/mol) of the lowest minima 
obtained for clusters composed of two pairs of PAH molecules. 
Circ. refers to the circumcoronene molecule. 



conformation 
type 


intertwinned 
stacks (a) 


single stack 
(b) 


parallel 
stacks (c) 


(coronene)4(HBC)6 


-1430.1 


-1500.7 


-1442.2 


(coronene)4(circ.)6 


-1725.7 


-1826.2 


-1749.9 


(HBC)4(circ.)6 


-2058.4 


-2176.2 


-2069.2 



TABLE IV: Binding energies (kJ/mol) of the lowest min- 
ima obtained for clusters of six large and four smaller PAH 
molecules, (a), (b) and (c) refer to the labels of Fig. 01 while 
circ. refers to circumcoronene. 



the two or four tips of the stacks. This result is again in 
agreement with the energetic arguments used to explain 
the trends in Table ITlTl 

In Table im we compare the relative energies of these 
three conformations for the three mixtures. Substituting 
the six innermost molecules of the single stack (b) always 
yields the optimal structure. This result reflects the extra 
stability of single stacks of large molecules, as seen pre- 
viously. Hence the geometry is essentially determined by 
the larger PAH molecule. This observation cannot easily 
be generalized to arbitrary numbers or small and large 
molecules, and it would be interesting to look further 
into the effects of changing composition on the relative 
stability of heterogeneous assemblies. 



IV. INTERMOLECULAR VIBRATIONS 

Our eventual interest lies in interpreting experimental 
infrared spectra by identifying the possible PAH assem- 
blies responsible. Although the present paper is mainly 
devoted to the prediction of structure, we wish to address 
here the intermolecular vibrations of the smallest clusters 
(dimers), leaving a more general discussion about larger 
assemblies and the couplings with intramolecular motion 
to future work. 

The normal mode frequencies were obtained by double 
numerical differentiation of the potential energy and ap- 
propriate mass weighting. The main vibrational modes 
are characterized in Table for all dimers of the same 
PAH ranging from pyrene to circumcoronene, as well as 
the zwittcrionic form of (coronene)2 for comparison. As 
expected from the strong anisotropy of the molecular in- 
teractions, the vibrational modes can generally be gath- 
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(a) 


(b) 


(c) 


(d) 


(pyrene) 2 


55.1 


5.5 


55.5 


3.6 


(coronene)2 


58.9 


4.3 


55.9 


5.6 


(ovalene)2 


59.9 


6.6 


69.6 


5.0 


(HBC)2 


62.1 


4.8 


58.1 


3.5 


(OBC)2 


62.1 


4.9 


64.5 


2.3 


(circ.)2 


62.9 


4.1 


61.0 


3.1 


coronene^ coronene ^ 


69.0 


6.8 


81.3 


47.4 



TABLE V: Frequencies (cm~^) for the vibrational modes of 
several PAH dimers. For the shearing (b) and bending (c) 
modes, an average of the two values is given. 



ered into two separate groups, the bending and breathing 
motions being the stiffest and the shearing and twisting 
motions the softest. The frequency of the breathing mode 
grows quite regularly with increasing PAH mass. For the 
circumcoronene dimer, this frequency is quite close to 
the value generally accepted for graphene layers, namely 
/o = 63.7 cm~^. This result indicates that our atomistic 
modelling is physically correct. 

The increase of the frequencies for larger PAH's is a 
natural consequence of the increasing strength of the 
repulsion-dispersion interactions. Both the repulsive and 
attractive parts of the potential increase roughly linearly 
with the surface area of the PAH molecules, because they 
are both short-ranged. At first order, the force constant 
scales approximately linearly with the PAH mass, there- 
fore the associated frequency remains constant. How- 
ever, the attractive part of the potential has a shorter 
range than its repulsive part, so it increases somewhat 
faster with the mass. Thus the effective force constant 
increases slightly more rapidly than linearly with the 
mass, in agreement with the observed increase in the fre- 
quency. This effect is magnified by the decreasing role of 
the Coulomb forces, which become negligible in the limit 
of two graphene sheets. 

The bending modes generally show similar frequencies 
to the breathing mode, and also the same pattern. How- 
ever, less symmetrical molecules, such as ovalene or OBC, 
have non-degenerate bending modes, resulting in smaller 
average values. Both the twisting and shearing modes are 
comparatively very soft, reflecting the small corrugation 
of the molecular planes. 

Lastly, the vibrations of the zwittcrionic form of the 
coronene dimer, chosen as a typical example, display sig- 
niflcantly stiffer modes. In this system, the Coulomb in- 
teractions become dominant at long range, which affects 
the shearing, and especially the twisting modes. 

Larger assemblies of PAH molecules show a greater va- 
riety of modes. Stack vibrations have been considered by 
Seahra and DuleyS in the one-dimensional approxima- 
tion of longitudinal modes only. However, the results ob- 
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tained in the previous section show that multiple-stack 
conformations should become more stable as more and 
more PAH molecules are added. This result calls for a 
more complete study of vibrational properties, including 
the possible influence of intramolecular motion. 

V. DISCUSSION 

A. Difficulties with optimization 

In the original basin-hopping global minimization ap- 
proach followed here, all molecules are translated and 
rotated simultaneously from their original position be- 
fore a new local optimization is carried out. Subsequent 
moves are then performed starting from the current lo- 
cal minimum in a Markov chain. This approach was not 
very efficient for large PAH assemblies and/or extended 
individual molecules due to the frequent intersections in- 
duced by the large amplitude motions. The parallel tem- 
pering strategy, which employs smaller physical moves, 
avoids such problems but can be quite slow in converg- 
ing toward low energies. In addition the set of tempera- 
tures needs to be adjusted for the configurations to be 
communicated to the colder trajectories. In practice, 
we supplemented the parallel tempering MC simulations 
with periodic quenches at the trajectory corresponding 
to r = 20 K. Several of the most stable structures were 
found by seeding the optimizations with a single main 
stack and adding the required extra molecules randomly 
in space. 

The results obtained in the previous section and the 
occurrence of small stacks as the main structural motif of 
larger assemblies suggest several ways for improving the 
basin-hopping technique in order to deal more efhciently 
with such anisotropic molecules as PAH's: 

(i) Randomly select a single molecule; choose its ori- 
entation from among those of the other molecules; 
move the molecule randomly in space to a location 
where it does not suffer any overlap. 

(ii) Select several molecules stacked together and move 
then simultaneously as a block, so that they do not 
overlap with any remaining molecule. 

(iii) Follow collective, large amplitude motion with a 
rescaling of all distances, such that the minimum 
pair distance is at least the PAH size; then choose 
random orientations for all molecules. 

(iv) Perform local moves that add or remove a single 
molecule at the tip of a stack. 

The above ideas can be combined into a single Monte 
Carlo approach, with appropriate probabilities for each 
type of move. They are suited to the specific case of pro- 
late molecules, which pose more significant difficulties in 
numerical simulations than spherical or oblate species. 
Such moves obviously constitute a bias to optimization. 



but they should improve the performance of the algo- 
rithm by guiding it toward the multiple-stack funnels. 



B. Influence of the potential 

The choice of the intermolecular potential is crucial in 
predicting the correct physical and chemical properties 
of macromolecular assemblies. This is especially true for 
smaller molecules, such as benzene, where small changes 
in the repulsion-dispersion or electrostatic parameters 
can easily lead to the wrong crystal structureiSi 

We did not find significant changes in our results on 
taking the Lennard- Jones parameters from other sources, 
such as the OPLS model of Jorgensen and coworkersi^ 
However, the same is not true for the electrostatic mul- 
tipolar contribution. The NBO charges obtained from 
the same DFT calculation on the coronene molecule are 
significantly larger in magnitude than the EPF charges, 
especially on the hydrogens f Table Hljl . The correspond- 
ing Coulomb energy is therefore much larger, and be- 
comes comparable to the dispersion energy. Such charges 
yield very different global minima, as illustrated in Fig. 
for (coronene) 13. The disappearance of the handshake 
structure. Fig. reflects the much longer range char- 
acter of the intermolecular interaction, which somewhat 
shields the molecular anisotropy. Here it becomes nec- 
essary to clearly distinguish between the anisotropy of 
the molecules themselves (the so-called shape anisotropy) 
and the effective anisotropy of their intermolecular poten- 
tial. Increasing the magnitude of the Coulomb interac- 
tion results in a stronger short-range repulsion, which is 
first manifested in destabilization of the simple stack, ei- 
ther crossed or parallel-displaced, for the coronene dimer. 
At their new equilibrium distance, the relative orienta- 
tion between two molecules is mainly driven by the inter- 
action between quadrupoles. For larger PAH molecules, 
we expect the partial charges and quadrupole to become 
less and less important as the carbon frame gets more 
homogeneous. On the contrary, smaller PAH's should be 
even more sensitive to the electrostatic efi^ects. This is 
consistent with the equilibrium structures, which favour 
the T-shaped motif over parallel stacks. 

In larger assemblies, isotropic interactions favor an 
icosahedral growth scheme over stacking. In the struc- 
ture displayed in Fig. |7|d, the molecular centers of mass 
indeed form a distorted icosahedron, similar to what has 
been found in (benzene) 13 by several authors iS2i2Si22i2£ 
Such structures are significantly easier to locate than 
stacked isomers, because the orientational degrees of free- 
dom can be neglected in a first approach. Here the third 
trick suggested in the previous subsection for improving 
the basin- hopping algorithm should be efficient. 

The MuUiken charges responsible for the icosahedral 
minimum are probably much too high and they do not 
support stable stacks of coronene PAH's, which is in dis- 
agreement with most previous studies. Conversely, set- 
ting the partial charges to zero does not alter the results 
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(a) 




(b) 




FIG. 7: Putative global minima obtained for (coronene)i3 
based on the same potential energy function, Eq. Q but dif- 
ferent sets of rigid partial charges obtained from (a) EPF 
charges, and (b) NBO charges. 



very much, thus confirming that stacking is mostly in- 
duced by the intrinsic molecular anisotropy. 



VI. CONCLUSION 

The work reported in this paper was devoted to locat- 
ing the global minimum for assemblies of large polycyclic 
aromatic hydrocarbon molecules. To do this we first con- 
structed an empirical potential from available data for 
the repulsion-dispersion part, as well as from some care- 
fully analysed density functional calculations for the elec- 
trostatic contribution. These extra ab initio calculations 
were mainly motivated by the lack of published data on 
the specific interactions between large PAH molecules. In 
an attempt to rationalize the various charge distributions 
obtained for the neutral and ionic PAH's, we adjusted a 
simple fluctuating charges model based on four param- 
eters. This model, used for reoptimizing the charges of 



the PAH molecules in the environment of their own sta- 
ble assemblies, indicates that the fixed partial charges 
approximation is quantitatively correct, at least close to 
the equilibrium geometries. 

Using a combined set of global optimization tech- 
niques, we have shown that clusters of PAH molecules 
grow by forming stacks. This optimization process 
turned out to be difficult because of the large amplitude 
random motions involved. Therefore the putative global 
minima we have found might be bettered in future work, 
although we expect the trends to be correct. We sug- 
gested ways of improving optimization algorithms dealing 
with strong shape anisotropy, which could also be useful 
for the liquid crystal phases of discotic molecules iS2*SiiS 

Stack formation appears to be mainly driven by the 
PAH size, since clusters of benzene and, to some ex- 
tent, naphthalene or anthracene, tend to keep some of the 
polytetrahedral order known in small atomic clustersiSS 
However, the one-dimensional growth of the stack is quite 
limited, and more close-packed forms become favored af- 
ter a certain size, which increases with the PAH diameter. 
These new forms are based on smaller stacks, and larger 
clusters are also made of the small stack motif. Even- 
tually, parallel stacks with alternate orientations of the 
molecular plane appear, as the precursor to the herring- 
bone crystalline arrangement. We also found that the 
electronic structure calculations of the partial charges 
had to be performed carefully, since overestimation of 
the Coulomb energy wrongly leads to more isotropic in- 
teractions, which destabilizes the stacks. 

The formation and physico-chemical evolution of 
PAH's is a question of great interest in astrochemistry. 
Recently, Rapacioli et alm^ have revealed the chemical 
link that exists between PAH molecules and small car- 
bonaceous grains. These authors suggested that the 
grains might be PAH clusters. The results presented here 
can provide guidelines to understand the nature and evo- 
lution of such systems in astrophysical environments. We 
note that all these results are only based upon energetic 
considerations and that entropic criteria have not been 
taken into account. Growth could proceed by monomer 
addition, but also by the addition of small clusters, pre- 
sumably as short stacks. Therefore, the aggregation ki- 
netics could also play an important role in the structure 
of the final aggregate. As we have seen, thermodynam- 
ics should favor segregation, since smaller PAH's are ex- 
pected to be located in the outermost parts of the aggre- 
gate. Upon addition of one or more extra PAH's, the re- 
arrangements required to reach the true global minimum 
involve very high energy barriers, due to the anisotropy 
of the molecules and their layered structures. In the in- 
terior of molecular clouds the temperature is very low, 
which might be an argument against the relevance of 
calculations based on thermodynamics only. However, 
the large timescales involved in molecular clouds, whose 
lifetime is close to ~10^ years, could allow the cluster 
to explore many configurations, eventually reaching low- 
lying minima. Another opportunity for rearrangement 
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into these stable configurations comes from heating by 
UV-visible photons originating from stars at the border 
of the clouds. It is in such regions that Rapacioli et 
al."^^ have reported photophysical data for PAH clusters, 
through their mid-IR emission spectrum and photoevap- 
oration properties. In these heated regions, we predict 
that for mixed clusters the smallest PAH's are located at 
the periphery and are therefore photoevaporated first. If 
this is the case, the composition of both aggregated and 
free interstellar PAH's is expected to strongly depend on 
UV processing. 

In this paper we also briefly discussed the intermolec- 
ular vibrational modes. Although our results are very 
preliminary, they already give an idea of the frequen- 
cies expected for such modes. For the stiffest modes, 
values between 55.6 and 81.3 cm~^ are found. They 
could correspond to bands in the astronomical spectra 
and can be searched for in the data obtained from the 



Long Wavelength Spectrometer, which was on board the 
Infrared Space ObservatoryiSi One difficulty is to disen- 
tangle these bands from the numerous intense molecular 
lines at moderate resolution!^ Hopefully, the PAH clus- 
ter features will be revealed by future space missions in 
the far-IR to submillimeter range, in particular, the Her- 
schel Space Observatory will explore the 16-175 cm~^ 
range. 

In addition to the intermolecular modes obtained from 
the present model, the influence of the cluster structure 
on the intramolecular frequencies needs to be quantified, 
and compared to available astronomical spectrai^ In- 
vestigating the coupling between inter- and intramolecu- 
lar modes requires a more realistic atomistic description, 
beyond the rigid-body approximation. Extensions of 
the tight-binding Hamiltonian developed by the Parneix 
group2& are currently under development. 
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